function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    w = e.getGaussWeight(0);
    xi = e.getGaussPoint(0);
    jac = f.jacXAt(xi);
    detJac = det(jac);
    B = zeros(3,6);
    B = e.computeB(f,B,xi);
    C = elements(i).constitutiveMatrix;
    if elements(i).isCrack==1
        %compute Xi
        due = ue - elements(i).ue_old;
        dxi = elements(i).Kxixi\(elements(i).renh-...
            elements(i).Kxid*due);
        elements(i).Xi = elements(i).Xi_old+dxi;
        %compute algorithmic tangent
            %the algorithmic tangent was not changed
        %compute strain
        epsilon = B*ue + elements(i).Gc*elements(i).Xi;
        %compute stress
        elements(i).stress = C*epsilon;
        %compute TGamma and update Kxixi
        lGamma = elements(i).lGamma;
%         elements(i).Xi
        n = elements(i).crack_normal;
        u_mu = elements(i).Xi'*n; %wrong here
        
        %modify Gc if u_mu got negative value
%         if u_mu(1) < 0
%             elements(i).Gc = -elements(i).Gc;
%             elements(i).Kdxi = w*abs(detJac)*B'*C*elements(i).Gc;
%             elements(i).Kxixi = w*abs(detJac)*elements(i).Ge'*...
%                 C*elements(i).Gc;
%             n = elements(i).crack_normal;
%             Ctilde = [elements(i).H*n(1) 0;
%                 elements(i).H*n(2) 0];
%             elements(i).Kxixi = elements(i).Kxixi+...
%                 lGamma*Ctilde;
%             disp(['Gc changed for element ' num2str(i)]);
%         end
        
        t = elements(i).ft+elements(i).H*u_mu(1);
        %TGamma = min(max(t,0),elements(i).ft)*n;
        TGamma = max(t,0)*n;
        if t < 0
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac);
        else
%             Ctilde = elements(i).H*[n(1)^2 n(2)*n(1);...
%                 n(1)*n(2) n(2)^2];
            Ctilde = elements(i).H*[n(1) 0;...
                n(2) 0];
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac)+...
                lGamma*Ctilde;
        end
        
        %compute renh
        renh = -lGamma*TGamma;
        renh = renh - w*elements(i).Ge'*...
            elements(i).stress*abs(detJac);
        elements(i).renh = renh;
    else
        %compute algorithmic tangent
        %the algorithmic tangent was not changed
        %compute stress
        elements(i).stress = C*B*ue;
    end
end
%
%
%
for i=1:ne
    e = model.getElement(i-1);
    face = e.getFace();
    if elements(i).isCrack==1
        %save Xi to Xi_old
        elements(i).Xi_old = elements(i).Xi;
    end
    %save ue to ue_old
    elements(i).ue_old = face.getUHatLocal();
end
end
